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ABSTRACT 


The problem of separated flow about a circular cylinder is studied 
analytically through the use of two potential flow models. Following 
a detailed review of the present state of knowledge, mathematical 
formulation of the models and appropriate computer programs are 
presented. These models are shown to be capable of simulating the 


asymmetric vortex shedding. 
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ENGLISH 


TABLES OF SYMBOLS AND ABBREVIATIONS 


Drag coefficient; defined by equation (1) 

Pressure coefficient; defined by equation (27) 

Radius of circular cylinder or sphere 

Differential operator 

Cartesian tensor notation for Eulerian total derivative 
Drag force on circular cylinder 

Lift force on circular cylinder 

Cartesian tensor notation for body forces 


Acceleration of gravity 
NEI 
A free parameter in Free Streamline Theory 


Mach number or number of free vortices 


Total number of free and attached vortices in a 
flow field 


Pressure 


Cartesian tensor notation for the three axial velocities; 





q = u; а, = у, апа q, = w 
Ц 
Reynolds number; for circular cylinder, Ве = т 
A 
Strouhal number; for circular cylinder, S - HT 


Strouhal period; U. T = length of wake generated in 
one period 


Time 
Velocity of the ambient flow; could be time dependent 


Instantaneous velocity in x-direction of a two or three 
dimensional cartesian system 


W(z) 


GREEK 


pe x 
УМ, 
~_ 


( 


NUM 


ЕЕ Б 


-- 
БЕ 


SUBSCRIPTS 


Instantaneous total velocity vector 


Instantaneous velocity in y-direction of a two or 
three dimensional cartesian system 


A function of complex variables; defined by, 


wlz)= Р(х,у) + P(x,y) 


Instantaneous velocity in z-direction of a three 
dimensional cartesian system 


Complex variable; z = x + iy, an undefined location; for 
example, z = x. T iy, could be the location of the nth 
vortex 


Strength or circulation of a vortex; defined as the líne 
integral of the tangential velocity component around a 
closed curve, or defined as the integral of vorticity 
over a surface bounded by a closed curve 


Specific weight; defined by, ἔ-ρ9 
Cartesian tensor permutation symbol 


Cartesian tensor vorticity vector; defined by equation 
(20) in Section II 


An angle measured from the positive x-axis 
Absolute viscosity 

Kinematic viscosity; defined by, 2% = + 
Density I 


Velocity potential function 


Stream function 


Cartesian tensor index; i 1,2,3; or indicates 
quantity pertains to an image vortex 


Cartesian tensor index; j = 1,2,3 
Cartesian tensor index; k = 1,2,3 
Vortex counter index; m = 1,2,.....,M 


Vortex counter index: n = 1,2,..... ,N 
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I. INTRODUCTION 


The postulation of an analytic solution for a viscous fluid flow 
about a circular cylinder has perplexed hydrodynamic investigators for 
several centuries. Although notes left by Leonardo da Vinci (1452-1519) 
indicate more than a qualitative interest, it was not until the early 
seventeenth century, during the evolution of Potential Theory, that the 
first formal attempts were made at quantitative analysis. Though there 
was some understanding of viscous effects during this period, they were 
basically qualitative and were not taken into consideration when in 1750 
D'Alembert found that: 


The force on a body in steady motion relative to an other- 
wise undisturbed ideal fluid is zero. [1]! 


This result, classically known as D'Alembert's Paradox, was of course 
unsatisfactory. Decades have elapsed since investigators found the 
qualitative reasons for the erroneous results given in this and other 
cases by Potential Theory. 

Philosophically speaking, to the engineer the utlimate goal of any 
analytic procedure is to provide a power of prediction for design ap- 
plications. As yet this capability is generally only realizable from 
experimental investigations and then only for certain flow characteris- 
tics. For a circular cylinder immersed in a steady cross-flow as shown 


in Figure 1, the drag force can be calculated from: 


И т”, 
με... 


λα... 2 
papel, ee (1) 


Lhe bracketed numbers refer to the item number in the List of 
References. 


ES 
ΝΡ 


Figure l. Circular Cylinder Immersed ia a Steady Fluid Cross- 
Flow and Indicating Relatienship to a Complex 
Coordinate System. 


C, is determined experimentally as a function of the Reynolds number as 
shown in Figure 2. Except for Re < 1, there is presently no general 
analytical method available to generate these results. This exception 
will be discussed in more detail in Section II. The primary point in- 
volved here is that for each geometric shape, or group of geometric 
shapes immersed in a steady flow field, there would have to be an ex- 
perimentally found C, curve for each case. In fact, for an incompres- 


D 


sible fluid, the form of equation (1) would have to be more general: 


rj 
ῄ 


f(Re), (2) 


and for a compressible fluid: 


rrj 
II 


Е(Ве, М), (3) 


For all but the simplest geometries, the definition of Re and M are in 
themselves difficult problems; not as to their general meanings, but as 
to specific combination of the physical parameters for a given applica- 
Exon. 

Though used only as an illustrative example above, the analytical 
determination of Ср has more subtle implications than those previously 
indicated. It requires a complete understanding of separated flow, 
particularly from the aspect of vorticity generation in the boundary 
layers befcre separation. It is this point that presents the greatest 
problem in the postulation of a mathematical model; for it is assumed 
that if a general analytical solution can be found for Ср in time- 
dependent flows, then analytic determination of all or most of the other 


flow characteristics will follow. Also, since a circular cylinder and 


its surrounding flow field can be easily represented in two dimensions 


11 
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by a function of complex variables, it is possible to easily transform 
the cylinder into other two-dimensional shapes using the technique of 
conformal mapping. For inviscid, incompressible fluids, this technique 
has been used to great advantage particularly in aerodynamics, however 
external knowledge of boundary layer effects was needed. 

So far, the discussion has centered around total drag. For future 
analysis it is necessary to break this into its two classical components, 
skin friction and form drag. The skin friction is the integral of the 
shearing stresses over the surface of the body and is dependent on the 
action of the fluid in the boundary layer. The form drag is equal to 
the integral of the forces normal to the surface of the body and is 
assumed to be dependent only on the inviscid fluid action outside the 
boundary layer. This assumption is made since an order of magnitude 
analysis of the Navier-Stokes equations indicates that the pressure 
gradient through the boundary layer is negligible as compared to other 
effects. The practice of considering the fluid outside the boundary 
layer as inviscid is due again to an order of magnitude analysis of the 
Navier-Stokes equations. This time, when considering only the region 
outside the boundary layer, the equations reduce to the inviscid equa- 
tions of motion and permit the application of Potential Theory to this 
region. Note that large Reynolds numbers have been assumed in both 
instances. 

Following the above reasoning further, leads to the conclusion that 
if no separation of the boundary layer occurs, the form drag is essen- 
tially zero since there is no disturbance to the external inviscid flow. 
At very low Reynolds numbers, this conclusion is reasonable and can be 
considered a basis for the famous "creeping motion" solution of G. G. 


Stokes in 1851. In this case skin friction predominates. At larger 
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Reynolds numbers, however, there is boundary layer separation and form 
drag becomes important. The external flow field is disturbed by regions 
of intense vorticity which causes a decrease in pressure behind the body 
as compared to the undisturbed inviscid solution. It is evident then, 
that without considering the direct shear force effects of viscosity as 
represented by skin friction, and the indirect effects represented by 
vorticity regions, the qualitative concepts reduce such that D'Alembert's 
quantitative solution results. 

Earlier in this discussion, mention was made concerning the practice 
of considering the undisturbed fluid outside the boundary layer as in- 
viscid. This requires amplification and a word of caution. In 1908, 

H. Blasius presented a solution to the boundary layer on a circular 
cylinder. As far as his paper was concerned the method of solution, 
rather than the result, was the most important contribution. Briefly, 

it involves the use of a series expansion for the free stream velocity 
just external to the boundary tavern If this velocity is taken as that 
found for incompressible, inviscid flows derived from Potential Theory 

on the surface of the cylinder, the angular distance between the upstream 
stagnation point, and the boundary layer separation point is found to be 
108.8 degrees; the center of the cylinder is taken as the Orrin ТЕ 
experimental velocity determinations are used, the solution yields a 


separation point at about 82 degrees. The observed separation point is 


For a complete presentation of this method, see [2], pages 154-162. 


The boundary layer separation point is defined as: oz gue 
where u is the fluid velocity parallel to the body surface, i 
and y is directed normal to the surface, with its origin at the surface. 


14 


at about 81 degrees. Herein then lies the problem with the above prac- 
tice: even at high Reynolds numbers there is sufficient interaction 
between the boundary layer and the external “inviscid” flow to cause 
considerable error in the "A ce of essentially exact methods of solu- 
tion. As will be shown in Section III, the separation point is critical 
in any method that attempts to mathematically generate a reasonably cor- 
rect drag coefficient curve. 

In recent years, extensive experimental work with the flow past 
bluff bodies has provided researchers with a clearer understanding of 
the mechanism of boundary layer separation. As a result, several models 
have been devised to portray the overall fluid action. However, to date 
none has provided the sought-after general Sol tionc. A survey of the 
literature indicates that there are two basically different analytical 
approaches being made. The first, and probably the most obvious, is the 
solution of the nonlinear Navier-Stokes equations using finite differ- 
ence techniques. The second is a potential flow model in which the shear 
layers are represented as point concentrated vorticies. Both these 
models, particularly the latter, will be described in greater detail in 
succeeding sections. 

Since this thesis does not propose a general solution, it is hoped 
that the following sections will help to clarify the problem and to 
provide a basis from which future work can be done. Section II reviews 
the literature on separated flows; concentrating on those aspects that 
relate to incompressible flow about circular cylinders. Section III 
presents the two potential flow models that the author used in an at- 


tempt to generate the aforementioned drag curves. 


4 
The development of the high speed, large storage capacity digital 
computer has also been extremely important to the study. 
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11, REVIEW OF THE PRESENT STATE OF KNOWLEDGE 


Because of the large number of references available on separated 
flows, the author has elected to limit this review through the use of 
the following guidelines. First, since there still exists a lack of 
basic understanding of the mechanism of separation only references to 
studies involving essentially incompressible fluid will be used. This 
enables the added complications of compressible phenomenon, both qualita- 
tive and mathematical, to be ignored. Secondly, a study of the literature 
reveals that several different geometries have been used both in experi- 
mental and analvtical work; the results not always consistent nor related 
to each other. Therefore, the circular cylinder is taken as the basic 
peometry for this review. The use of a circular cylinder implies, to a 
large extent, two-dimensionality; which is another simplification. 
Though particularly helpful ΠῚ a qualitative standpoint, this simplifi- 
cation is also useful mathematically. Thirdly, the information will not 
necessarily be presented in chronological order as is so often done in 
reviews of this nature. It is felt that a loose division by Reynolds 
number is more pertinent since the analytical solutions thus far derived 
are only applicable to given ranges of this "magic number," as von 
Karman calls it. Lastly, the discussion will in general be limited to 
the laminar flow regime. As the Reynolds number increases above the 
critical value, increasing consideration must be given to compressibility 
and eventually to cavitation; two complicating subjects outside the scope 


Of this thesis. 


А. REYNOLDS NUMBER LESS THAN ONE 


In this Reynolds number range, the viscous forces acting on the 
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cylinder are considerably larger than the inertial forces. This leads 
to the conclusion that if the inertia terms are dropped from the Navier- 
Stokes equations, an approximate solution may result that would compare 
favorably with experimental data. Stokes, as previously indicated, did 
this in 1851 when he presented the first solution for the drag on any 
body. Though he used a sphere in his solution, it is still worthwhile 
presenting here. The mathematical details are not the same, but the 
basic steps are. 

The complete Navier-Stokes equations for incompressible, constant 
temperature conditions are (expressed in three-dimensional Cartesian 
tensor notation [4]): 


glas Mon BN . 


The incompressible continuity equation is: 
Fie =O (5) 


By dropping the inertia and body force terms, then cancelling density 


from both sides of the rearranged equation, one has 
= ^ o e 
PaM hasis © 


The four equations represented by (5) and (6) have the same usual 
boundary conditions used for any boundary layer problem of this nature 
when the body is immersed in an infinite fluid: the fluid velocities 
normal to and tangent to the body surface are equal to zero. 

After taking the divergence of both sides of equation (6), it 


becomes: 


iA dig ΄ (7) 
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But the right hand side of equation (7) can be written as l р » .. 
š ^^ ere 2224 

which is the second derivative of the continuity equation and there- 

fore equal to zero. This means that in creeping motion, the pressure 


field satisfies Laplace's equation and hence is a potential or harmonic 


function. Now, by introducing the two-dimensional stream function, 


eT ue j (8) 


into the right hand side of equation (7), dropping o and μέ the 
244. 


two-dimensional biharmonic equation results. 
Ys 2 d 
Dd DX) (9) 


Without going into the mathematical details, it is sufficient to 


say that the pressure and velocity components reduce to: 


2 отн а a 
бере ) (+= + | 


(10) 








The geometric configuration is such that the center of the sphere is at 


18 


the origin of the coordinate system and is placed in a positive velocity 
field, Uo , which is parallel to the x-axis. 

Next, by using the standard definitions of shear stress along with 
equation (10), the viscous and pressure forces can be found by integra- 
tion over the sphere's surface. The total drag for the sphere then 


becomes, 
E-6TAcCH, ° (11) 


If a drag coefficient is formed in the normal manner (i.e., referring 


to the dynamic head and frontal area), then, 
2 
2/1 
= Wc (22 Ц ) e 12) 
Equating (11) and (12), then, 


С = 24 (13) 
А Re 
Comparison of this result with experimental data indicates that 
these equations are valid only for cases where Re < 1. A classic case 
of application is that of small, light particles, assumed spherical in 


shape, settling in a fluid. Assuming that the particle is at its ter- 


minal velocity, then summing the forces yields, 


2 2 i 
9g С E Ж шз E Ц. z 
Letting â Y= εί oe шд) then equation (14) reduces to 
εδ. 9 


= — 15 


which is the famous Stokes' number. 
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C. W. Oseen, in a paper published in 1910, improved on the above by 
taking the inertia terms of the Navier-Stokes equations partly into ac- 


count [2]. By doing so, the resultant drag coefficient became: 


е (1+2 Re) 2 (16) 


Comparison of (16) with experimental data shows favorable correlation 
up to a Reynolds number of about 5. 

In recent years, other investigators have formulated other solutions, 
however, their main contributions have not necessarily been to engine- 
ering but more to mathematical prowess.  Kaplun [5] based his on an 
asymptotic solution of a series expansion of the Navier-Stokes equations. 
The results are generally restricted to two-dimensional bodies with 
equivalent Reynolds numbers less than one. Sir Horace Lamb [6] also 
presented a solution applicable to this range. He followed the basic 
ideas of Oseen, however, where Oseen used a sphere, Lamb used a circular 


cylinder. 


B. REYNOLDS NUMBER BETWEEN THREE AND FORTY 

In this Reynolds number region, the boundary layer is still laminar, 
however, separation occurs causing the formation of two vortex cores in the 
wake of the cylinder as indicated in Figure 3. This region is sometimes 
referred to as the 'Twin-Vortex Stage.' А great deal of experimental evi- 
dence indicates that the cores are generally stable in this region pro- 
vided that no external excitation occurs. It is in this region also 
that the greatest successes have occurred with the "exact" solution of 
the Navier-Stokes equations. Remembering that, above a Reynolds number of 


about one, the inertia forces are too large to neglect, the phenomenon 
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requires the solution of the full Navier-Stokes ООО оз.” This is most 
readily done through the use of either the forward or backward finite 
difference techniques. Until recently, without the advent of the digital 
computer, these methods have been difficult to employ. Thom [7] in 1933 
used a relaxation method of hand calculation to produce the streamline 
pattern shown in Figure 4 at a Reynolds number of 20. 

It is true that the relaxation method can be programmed for a digi- 
tal computer and made to converge to any realistic degree of accuracy, 
however, it is also true that computation time increases rapidly for a 
slight increase in desired accuracy. Thus this method is not considered 
the most desirable. The Gauss-Siedel iteration method is somewhat 
different in procedure, but can require the same amount of computer 
time for solution. [It should be noted that both of these methods can 
be modified if done by hand since one can make shortcuts that are diffi- 
cult to program for the computer. Even so, the digital computer reduces 
the time required from ".....about one year and a half with twenty work- 
ing hours every week, with a considerable amount of labour and endurance 
o cos "[8] to a total computation time of about three to five hours. 

Even this could probably be reduced by streamlining the iterative proce- 
dures. Robertson [l, р. 366] summarizes these methods quite well. 

Though these analytical methods produce interesting flow patterns, 
their primary reason for use is to give needed information about sepa- 
rated flow. Unfortunately, they do not appear to explain the mechanism 


of separation or more importantly the feedback aspects of the wake on 


5 

Body forces, f., are usually neglected in these solutions since 
their presence adds nothing to the understanding of the problem and in 
general add unwanted complications to the mathematics. 


2 2 
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the vortex generation in the boundary layer. The experimental and ana- 
lytical work does yield approximate Reynolds numbers for the beginnings 
of separation and the formation of the Kármán vortex street; about 5 and 
40 respectively. From an analytical standpoint, the mesh size of the 
finite difference grid and the method of defining the conditions along 
the outside grid DOR. can affect the resultant separation and wake. 
Experimentally, just the placement of the confining walls of the ap- 
paratus used for testing can increase the Reynolds number, at which the 
wake becomes unstable, by more than 300 percent [9]. It should be noted 
here that although the wake is generally considered to be stable for 
steady flows in the Reynolds number range under discussion, it can be 
made unstable both theoretically and experimentally under artificial 
excitation. The ease with which this can be done, of course, increases 
as Re approaches the range of instability, making exact determination 

of the Reynolds number at which the vortex shedding begins very diffi- 
cult. This same idea holds true for the determination of the separation 
Reynolds number. In both cases there appears to be some experimental 
evidence to indicate that the Reynolds number is not the exact similarity 
parameter that should be used in comparing results; the reasons for this 
are not entirely known. It is the belief of the author, however, that 
these disparities in results are probably due to variance in technique 
rather than to any fault of the parameter. It is recognized that in 
actuality there may be a three-dimensional effect that has not been taken 
into consideration theoretically or experimentally. Herein may be much of 


the reason for the differences in experimental results. 


6 

Requiring that the velocity on the surface of the cylinder be 
zero provides the boundary condition for the grid points on the 
cylinder surface. 
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Figure 5. Beginnings of Vortex Street, Ве ^ 40. 


Figure 6. Karman Vertex Street Indicating Cyclic 
Shedding Vertex Cores, 80 < Re < 300. 
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E REYNOLDS NUMBER BETWEEN FORTY AND THREE HUNDRED 

The primary phenomenon represented in this region is the Karman 
vortex street; from the condition of an almost stable wake with slow 
viscous decay to the ultimate condition of cyclic detachment of vortex 


cores. Figures 5 and 6 indicate the flow patterns at these two "ex- 


tremes." 


This region is also characterized by a smooth decrease in drag 
coefficient and an almost linear increase of Strouhal number, S. 

Though study in this region has not produced any startling insight 
into the mechanism of separation, it has provided much information about 
vortex streets. Basically, it has been found that given two parallel 
vortex sheets separated by some initial distance, a solid body need not 
be present for a vortex street to develop [10]. The creation of these 
sheets is of course dependent on the viscous boundary layer effects 
caused by the presence of the solid body. 

Due to the oscillatory nature of the vortex core separation into 
the wake, there is a cyclic shift in the separation points of the 
boundary layers on either side of the cylinder. This action is evidence 
of the strong "feedback" effect that the wake conditions have on the 
boundary layer in this Reynolds number region. This condition exists 
during the so-called steady state period after the various parameters 
such as drag coefficient, Strouhal number, etc., have settled to fairly 
constant values. In impulsively started flow situations, the flow pat- 
terns above and below the cylinder are basically mirror images of each 
other. The mechanism that causes eventual asymmetry is not clearly 
understood. In recent years, with the help of the digital computer, 
this asymmetry has been generated in finite difference solutions by the 
artificial insemination of small amounts of extra vorticity into one of 


the boundary layer separation regions, once the twin core pattern was 
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fairlv well developed. Using this method, Fromm and Harlow [11] were 
able to duplicate the flow patterns around rectangular cylinders at 
these intermediate Reynolds numbers. At the lower numbers the system 
"absorbed" the extra vorticity with only minor disturbances to the wake, 
while at higher Reynolds numbers, the vortex cores separated cyclicly to 
form a vortex street. 

Though the details of solution of the Navier-Stokes equations varies 
with the investigator, the same general procedures are used by all. The 
following is an attempt to clarify the salient points involved, but 
first, the reason for using these equations should be reiterated. 
Through this Reynolds number region, there appears to be a balance of 
sorts between the viscous and inertia forces, or effects, which implies 
that neither can be neglected or ignored. In the potential model that 
will be discussed later, the viscous forces are neglected from a drag 
standpoint. In the Stokes solution for Re € l, the inertia forces were 
neglected with no apparent loss of practical accuracy. This region then 
appears to be nature's transition region from one extreme to the next, 
and therefore the complete governing equations must be used. 

The method, in general, is as follows; dropping the body force 
term from equation (4), and expanding to the two-dimensional case, 


there results: 


ο 


LL 1 ӘР _ 205 Xu 
+ V — + ru oS y ( 3X» | (17) 


ajay 
+ [> 
Q. 
Viv 
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and, 


πι 
|< 
4- 
E 
kf 
+ 
< 
9; 
< 


Eliminating pressure by cross-differentiation and subtraction, one has: 


(19) 


The vorticity vector is defined as 


c c | 20 
“Z, f< ¿ co 
which means that since only two dimensions are involved, vorticity 


exists only parallel to the Z-axis. Therefore, letting 
oV p 
C=C = мо (21) 
3 OX p 


and substituting back into equation (19), one has 


ии D 5 XE 
ΕΜ μα μι... me Oye ) 


or finally, 





7 Р 
u, V, and P must be assumed continuous functions of x, y, and t. 


The system of equations needed becomes complete by introducing 


the stream function relationship, 


T. m ^a hj (8) Repeat 


into the continuity equation, LL = 0. For the irrotational case, 
this would result with the Laplacian of the stream function equal to 
zero. However, there is rotation present and this introduction must be 
handled differently. By adding 2 to both sídes of the continuity 


equation, one obtains 


ди оа а дмс мые 
bm 8) Mor P 5 (23) 


The rearranging of the left side produces; 


д д ВРИЕ 2 
ЕКА 2 (24) 


Now substituting for u and v the appropriate expanded relation from 


av 
equation (8), i.e., u = 9 and v = - —— , and then performing 
044 dX 


the partial differentiation as indicated in front of each bracket, the 


equation becomes: 


W. СС? 


This of course reduces to a form of Poisson's equation, 


ee 27679 (25) 
4х? dy” 


Equation (22), when written in terms of V and & reduces to 


оф 9х _ д ох $i. m 
Jt “ΣΣ dx Әх dogs = >| ° Lx) ow 
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which is somewhat easier to program since the necessity of finding the 
velocity gradients has been eliminated. If desired, these could be 
found from the other forms of the equations given above. As before, 

the assumption concerning continuous functions noted in footnote 7 still 
applies. 

The specific finite difference formulas used to convert equations 
(25) and (26) into a numerically solvable system depends basically on 
the largest allowable truncation error and on the mesh size of the grid 
points. Crandall [12] discusses quite thoroughly all the aspects in- 
volved for both these considerations. It is sufficient to say that 
greater accuracy requires less truncation error and a finer mesh size. 
Generally, since finer mesh sizes require larger computer storage cap- 
abilities, they are used only in the areas of the flow field where large 
gradients are expected. 

There are in general two basic methods of solving this system of 
equations. The one most widely used, until recently, has been of an 
iterative nature such as the relaxation or Gauss-Siedel methods; or 
modifications thereto. The primary concern here is that of convergence 
and stability of solution. The Gauss-Siedel method is adversely sus- 
ceptible to both of these while the relaxation method generally does 
not have a stability problem, but can require an abnormally large number 
of iterations before some degree of convergence is attained. These 
methods are however, fairly simple to program. With the advent of the 
digital computer, more sophisticated numerical methods can be used to 
solve the system of equations explicitly using matrices. Methods such 
as Gaussian elimination with back substitution can produce a very ac- 


curate inverse for a matrix that is not singular, and a properly 
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constituted matrix representing a physical system should not be singular. 
For large systems of equations, the coefficient matrix can be reduced to a 
triangular form through the use of a sequence of linear operations. The 
inverse of this reduced matrix is then used rather than the inverse of 

the complete system. 

Besides Fromm and Harlow [11] several other investigators have 
studied this Reynolds number region using modifications of the above 
procedures. Most numerical solutions have involved rectangular cylinders 
due to the ease with which Cartesian coordinate finite difference equa- 
tions can be applied. Thompson [13] used a Gauss-Seidel method but 
appears to have assured convergence through the use of an “acceleration 
parameter" which is determined from a consideration of the eigenvalues 
of the matrices involved. His results compare favorably with experi- 
mental data, however, some of his computer runs required as much as nine 
hours on an IBM 360/40. This could have serious implications as far as 
accuracy is concerned, since each time step must adequately converge be- 
fore the conditions for the next time step can be calculated. If con- 
vergence presents no problem, then there is another important result of 
his work that should be mentioned: no artificial introduction of vorti- 
city is reauired to cause the eventual formation of the vortex street. 

Wang [1^] used cylindrical polar coordinates in a matched asympto- 
tic expansion solution of the Navier-Stokes equations. Though the 
results offer no new insight, the method is indicative of the mathe- 
matical attempts being made.  Ingham's paper [15] offers constructive 
criticism of Payne's method of finite difference solution [16], and is 
indicative of the pitfalls waiting for the unwary. Ingham's examples 
of variation of calculated C, for different mesh sizes at the same 


D 


Reynolds numer is quite dramatic. 
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Q, Degrees measured from forward stagnation point. 


Figure 8. 


Experimental and Computed Pressure 
Coefficient and Computed Surface 
Vorticity (5 ) at Low Reynolds Numbers. 


A Computed Ср» Аре! с, Re = 40 

B Computed Cp, Kawaguti, Re = 40 

C Computed 5 , Apelt, Re = 40 

D Experimental, Linke, Re = 2800 

E Computed Cp, Potential Flow Theory 
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The pertinent experimental and analytical data derived by investi- 


gators in this Reynolds number range are shown in Figures 7 and 8. 
Important contributors of this information include Taneda [17], Tritton 
[18], Roshko [19], Rosenhead [20], Wille [21], and Apelt [22]; besides 


those previously mentioned. 


D. REYNOLDS NUMBER BETWEEN THREE HUNDRED AND CRITICAL 

This subcritícal region is in general characterized by a relatively 
constant mean drag coefficient and a progressive loss in the distinct 
visualization of the Karman street. As the Reynolds number increases 
there appears to be less "feedback" effect of the wake on the separating 
boundary layers which are still laminar. The Strouhal number has a 
tendency to decrease slightly with an increase of Re. It is in this 
region that the potential flow point vortex model of the feeding sheet 
appears to have the best possibility of providing information since the 
inertia forces are so much larger than the viscous forces. The use of 
the finite difference Navier-Stokes equations to provide a solution 
seems to break down due to the increased turbulence in the wake. Al- 
though a much finer mesh size might alleviate this problem, the fact 
that the equations now being used do not provide for velocity perturba- 
tion is probably the more basic reason. There appears to be no argu- 
ment that the equations are applicable for use in describing the boundary 
layer on the front part of the cylinder almost up to the critical Re of 
about 5 x ο As the Reynolds number increases from the low end of 
this region, the boundary layer effects appear to feed vorticity into 
the wake at an ever increasing rate causing the wake to become extremely 
turbulent. As this occurs, the well ordered vortex core regions of the 


Karman street essentially disappear. Although the Strouhal number can 
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still be measured, the dominant periodicity seen in the lower Reynolds 
number wake flows all but disappears. Experimental studies using turbu- 
lent flow measuring techniques have confirmed the increasing turbulence 
with increasing Reynolds number. Ironically, the feedback problem di- 
minishes, but the problem of transition from laminar to turbulent flow 
conditions becomes more important. The understanding of the latter 
problem is even less than that of the former. The end result is that 
the applicability of the laminar Navier-Stokes equations decreases with 
increasing Re. 

In the "critical" Reynolds number region, loosely defined as 2 x 10? 
< Re < 109. there is a sharp decrease in the mean Cy down to about 
0.3.° According to Schlichting [2, p. 39], this is due to the following 
reason: when the total drag on any body is primarily form drag, the 
transition from a laminar to a turbulent boundary layer causes a decrease 
in the size of the wake and results in a decrease in drag. This can be 


interpreted as saying that the drop in C, is a boundary layer phenomenon, 


D 
though it does not quite explain the entire mechanism. After this transi- 


tion occurs, C, increases to a value slightly lower than that before 


D 
transition. There is also a return of a distinct vortex street. This 
time however, the ''vortex cores" are highly turbulent internally. As 
the Reynolds number increases beyond this region, compressibility ef- 
fects become greater, further complicating the problem. Depending on 


the fluid, very high Reynolds numbers can bring on the added, sometimes 


predominant, problem of cavitation. 


Ora M 

This is true for a hydraulically smooth cylinder. The roughness 
of the cylinder surface can have a major effect on this minimum mean 
drag and the Reynolds number at which transition occurs. See [2, 
2.0621]: 
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О, e-rees measured from forward stagnation point. 


Figure 9. Experimental Preesure Coefficients Found 
at Medium and High Reynolds Numbers 
Indicating Laminar (S1) and Turbulent 
(57) Separation Peints. 


A Linke, Re = 2.8х103 D Flachabart, Re = 6.7х105 
BíFage and) Re = 1. 1x10? E Roshko, Re = 8.4x10° 
clFalkner Í Re = 2.1x10? 


F Potential Flow Theorv 
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Experimental studies of the pressure distribution on the cylinder's 
surface have yielded some valuable information necessary for future 


analytical postulation. The pressure coefficient, defined as 


TT P. 27 
C. Sre (27) 


is an excellent parameter for theoretical and experimental comparison 
particularly for creeping motions (Re < 1), and at Reynolds numbers where 
the boundary layer separation points maintain fairly constant positions. 
At very low Reynolds numbers, the comparison is quite good especially on 
forward part of the cylinder, between points A and B in Figure 1. Be- 
tween points B and C, the actual pressure decreases somewhat since energy 
is expended in the boundary layer; i.e., not all the kinetic energy 
acquired from À to B is recovered as potential energy as the fluid de- 
celerates toward C. In Potential Flow Theory, points A and C are of 
course stagnation points, and whereas the pressure at À can be well ap- 
proximated by this theory over the entire Reynolds number range, the 
pressure at point C cannot as indicated in Figure 9. Note how curve D, 
found for a Reynolds number above the critical, more nearly approxi- 
mates curve F (potential flow). Photographs of the wake pattern pro- 
duced for this post-critical Reynolds number shows a narrow wake with 
the rest of the flow field similar in appearance to that of an ideal 
fluid. Note also that after the boundary layer has separated, оң as- 
sumes a nearly constant value, not the same for all Reynolds numbers 
however. 

Though ínformation about the Strouhal number is still incomplete 
above Re ο 102, there is an interesting similarity between the inverse 
of this parameter and C.. The curves for the two values, plotted over 


D 
the whole range of Re from the beginning of core shedding through the 
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critical region, are very similar in shape. In fact, the similarity is 
so good that the drag coefficient could be written as 
1 
EE s^ 1.4 

with the resulting value having an accuracy within about 15 percent. 
The explanation for this result is not complete, however, the point here 
is that there is apparently a strong conceptual connection between core 
shedding frequency and drag. It is interesting to also note that through 
the region of nearly constant drag, not only is 1/S fairly constant, but 
that the oscillations of the boundary layer separation points are gen- 
erally quite small in amplitude. 

Investigators have attempted several approaches in the search for 
an analytic solution in this region and there have been some notable 
successes. Each approach, however, involves the use of at least one 
experimentally determined or arbitrarily selected parameter. 

The free-streamline theory developed by Kirchoff, idealizes the 
free shear layers that separate from the cylinder by considering them to 
be surfaces of velocity discontinuity; hence a free streamline. This 
supposition works quite well for liquid jets discharging into a gas, 
however, for cylinder wakes, because of the assumption that the pressure 
in the wake is the same as the free stream pressure, the resultant cal- 
culated drag coefficients are much lower than those experimentally observed. 
The agreement can be greatly enhanced by the introduction of an arbitrar- 
ily selected parameter usually symbolized by K. This parameter is con- 


cerned with the depth of discontinuity (notch) in the logarithmic 


9 
hodograph plane, and is therefore concerned with the ratio between the 


το this plane, the two-dimensional physical flow is represented by 
straight lines parallel to the two coordinate axes, and is the result of 
conformal mapping. See [1, p. 351]. 
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free stream-line velocity and the free stream velocity. For very high 

Reynolds number ποσα where cavitation exists, K is fixed by the cavi- 
tation index, 0 . Fairly good agreement between experiment and theory 

exists in this case, particularly with respect to cavity flows. 

Dennis and Shimshoni [26] have probably been two of the more suc- 
cessful investigators in analytically deriving the mean drag coefficient 
over the entire laminar boundary layer Reynolds number range. Their solu- 
tion involves series RE À of the Navier-Stokes equations, set up 
for finite difference evaluations. At slow, creeping motions, their 
solution converges to that of Stokes and Lamb, but at Reynolds numbers 


above about 1.5, the values of C, are about 5 to 20 percent higher. 


D 
They allow no asymmetry of vortex shedding. The mesh size they use ap- 
pears to be much smaller than that of other attempts. 

As noted earlier in the discussion, this Reynolds number region is 
characterized by the predominance of inertia forces over viscous forces. 
Also noted was that this leads to the possibility of using point con- 
centrated, potential flow vortices to represent the shear layers or 
feeding sheets. Among those who have worked with this approach, 
Sarpkaya [35], Bryson [29],and Gerrard [25] have made some notable 
strides. The description of the work of the first two will be left to 
Section III, since it is basically with their ideas that the author 
worked. Gerrard's work, though similiar in mathematical approach to 


the others, appears to indicate a greater concern for wake characteris- 


tics than for finding solutions to the basic separation problem. 


и to Oseen and Lamb. 
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Besides the basic points covered here, there is a great wealth 
of experimental information on the more detailed aspects of the sub- 
ject that should be covered by any reader desiring to proceed further 
with the investigation. Toward this end, the author has attempted to 
compile the more notable of these in the List of References. Although 
some references have not been mentioned specifically herein thus far, 


this is not intended to mean they are of lesser value. 
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III. TWO ANALYTIC MODELS 


A. INTRODUCTION 

In the previous two sections, the discussion has centered on those 
few basic known facts and some of the many problems involved in this 
study. This section will present, in detail, two different analytic 
models, based on Potential Flow Theory, which are theoretically ap- 
plicable to the relatively higher Reynolds region discussed in Section 
II (D). 

The first model, called Model A hereafter, is conceptually more en- 
lightening than the second; hence the reason for its presentation first. 
It is basically derived from Bryson's [29] model, but where his model 
was limited to the growth of axisymmetric vortex cores, the equations for 
Model A are recast to allow "forced" asymmetry to be included into the 
solution. To the author's knowledge, this is the first time this has 
been done. This method, like almost all the others, has no provision in- 
cluded in the governing equations for "natural" asymmetry. Thus if asym- 
metry is to be studied, its beginnings must be forced. In brief terms 
then, the basis of this model is that large vortex cores are fed vorti- 
city through an infinitely narrow feeding sheet, or shear layer, con- 
necting the center of the core with assumed stagnation points on the 
cylinder's surface. 

In the second model, Model B, the shear layers are represented by 
a string of small strength vortices, called a vortex sheet; see [1, р. 
111]. The equations for this model were also derived and programmed to 
allow forced asymmetry. Whereas Model A represents a continuously fed 


large vortex core, this one creates numerous small vortices which are 
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allowed to roll-up, coalesce and/or cancel at will. There are four birth 
areas in this model, two for the primary vortices, and two for the sec- 
ondary СРР: This arrangement is a modification of Sarpkaya's 
model [35] in which there are only two axisymmetric birth areas.  Com- 


parison between the two will be left until later in this section. 


Dx BASIC POTENTIAL FLOW EQUATIONS 

The basic general equations used in both models will be presented 
in this subsection. The actual application will be left to the dis- 
cussion of that particular model. It is not felt that the history be- 
hind these equations is relevant to this thesis since the interested 
reader can refer to Sarpkaya [36] for this information. 

As mentioned earlier, in Section I, a two-dimensional flow system 
can be represented mathematically by equations written in complex vari- 
ables. Since the mathematical development of this representation is 
quite lengthy, the following is only a brief part of what can be found 
in many reference books and papers on hydrodynamics. To start with, 
it is well to review the basic assumptions. First, the fluid is con- 
sidered inviscid. This implies that there are no viscous forces pre- 
sent and therefore that the models derived should be applicable to the 
higher laminar Reynolds number range. Secondly, since this aforemen- 
tioned Reynolds number range is generally devoid of compressibility ef- 


fects, ~ will be assumed constant. The usual "steady flow" assumption 


11 
Vorticity is shed from the cylinder in the regions of boundary 


layer separation. These regions are associated with the "primary" 
boundary layers on the cylinder's upstream side and the "secondary" 
layers due to fluid back-flow on the downstream side. Since vorticity 
is created in these regions, they are called 'birth' areas, and the 
small strength vortices created are called "nascent" vortices. 
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will not be made, because all the flow characteristics appear to be 
dependent on the manner in which the boundary layer generates vorticity 
and on the manner in which it diffuses into the wake. The key problem 
then is the formulation of a method by which vorticity can be generated 
in a supposedly inviscid fluid. This is in direct violation of Kelvin's 
theorem for inviscid fluids which says that vorticity once created, can- 
not be destroyed, and if not already existing, cannot be created. In 
order to overcome this, boundary layer information must be taken into 
O aeration. 

From basic considerations of hydrodynamics, an undefined number of 
vortices outside the boundaries of a circular cylinder, suspended in a 


uniform flow field as in Figure 10, can be represented by 


Wiz) »-UL(z« )* | S E ix(-z) 
N P (28) 
e 
= 2 |, n ΓΝ | ë 


There are no symmetry restrictions on this equation. The normally ex- 
cluded terms have been left ο... Next it can be shown that for a 
general location Z , on or outside the cylinder, the fluid velocity is 


given by: 


А 
4 (2) у гу pt (S -)«g x 2 δ - 5 - Uu (29) 











dz "2 E л n=} Z 


2, 





12 
These are the terms that are excluded by virtue of the condition 


that the vortices shed from the cylinder leave a circulation opposite 
to their own on the cylinder. бее [38, p. 364]. 
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Figure 10. Circular Cylinder in the Presence 
of Numerous Vortices and a Uniform 
Flow Field. 
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Now that if 2-2 , the velocity goes to infinity. This is to be ex- 
pected since the velocity at the center of an ideal, inviscid vortex 
is equal to infinity. This same equation can, however, be used to 
find the absolute velocity of the vortex with respect to cylinder, here- 
after called the convective Or ter velocity, by letting 257, and dis- 
regarding the term that goes to infinity. This is due to the fact that 
the motion of a vortex is governed not by its own internal velocity dis- 
tribution, but by the forces induced upon it by other entities in the 
flow field. See [1, p. 124]. 

Next, the drag and lift forces on a body can be calculated from the 


generalized Blasius' therorem, 


: 4 qdW( 
ον 


Writing 


dz edo фи E S 








Б+2Е = (Х,-4У)+(Х,+23) о οὐ 


equation (30) can be broken down into its so called "steady state" 


2 


and the time dependent part, 


X, 4X, =ip 3p Wade Е (33) 


The steady part, using equation (29), can be evaluated [35] such that: 


N 
X tiy Sip S Dln tip) e en 
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The evaluation of equation (33) must wait until the discussion of the 
individual models since it involves some assumptions as far as the 
method by which E is to vary over the feeding zones. To continue, 
non-dimensionalizing can easily be done by referring all the length, 
velocity and time quantities to U and c; or for example, 
LSU E ) 
ЕЕЕ (35) 


Г-и. е 





апа, 
/ 
Е = < 
С... 
Substituting equations (35) into (29); letting U = 1.0; c = 1.0, 


and dropping the primes for the sake of simplicity, this becomes; 





N N 
с N DNI CNN Г 
ку = | | 2? — -$ - (36) 
d ^ mec N=) Е — 
y. 


Up to this point a discussion of image vortices has not been 
necessary, however, certain relations will be needed for future de- 
velopment of equation (33) in succeeding subsections. It can be shown 
by interpretation of Milne-Thompson's circle theorem (38, p. 157], or 
by using the method of images [38, p. 220], that the second summation 
in equation (29) and consequently in (36), are in reality "image" 
vortices inside the cylinder. The location of the image vortex, 

| „ЭИ у | | 
corresponding to the real vortex at 2 18 Z . Its circulation 
is in the opposite sense as compared to that of its corresponding real 
vortex. Positive numerical values for circulation indicate rotation in 
the counterclockwise direction, while negative numerical values indi- 


cate clockwise rotation. In the derivation of the equations, the 
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symbol Е has no inherent sign and can therefore take on positive or 
negative values as the situation warrants. 

Since there are conditions well described in [38] that make equa- 
tion (36) not applicable to the internal area of the cylinder, Sarpkaya 
[36] shows that the velocity of the image vortices are strictly a func- 
tion of the position of the image vortex and the corresponding real 
vortex's velocity. By manipulating his expressions, the following equa- 


tion for the image velocity results: 


5) V, = OL 





(37) 
zy 


To reiterate, the term whose denominator equals zero is disregarded 


when solving for sus Е іу 


C; MODEL A 

As previously mentioned, this model is primarily a modification to 
Bryson's solution [29]. Though his solution appears to give reason- 
able drag coefficients during the beginning time periods of an impul- 
sively started flow condition, once the peak is reached, the drag 
coefficient drops off to fluctuate about zero; Figure 11. This of course 
is not correct as Sarpkaya's superimposed experimental data [33] indi- 
cates. Though Bryson stops his solution when the vortices begin feeding 
circulation back to the о аан this in itself is not enough to 
produce a reasonable drag curve. It was hoped that by studying his 


method, modified to allow asymmetry, more could be learned about the 


e 
The condition occurs when m (to be developed later) for a 
given vortex changes sign. If the solution is not stopped at this 
point a negative drag coefficient results as shown in Figure ll. 


τ» 
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effect of vortex core shedding on the drag coefficient. Before a 
discussion of those results in undertaken, the equations must be 
developed. 

Using Figure 12 as the basis for nomenclature definition, the 
equation development starts with the necessary modifications to equa- 
tion (36). Though these equations can be constructed to allow more 
than two vortices, thís will not be done at this time in order to show 
a simplified version of the technique. As will be indicated later, the 
full equations are οῦ δε Неа by the simple addition of a summation term 
to each of the following equations. So, based on two vortices, the 


velocity of vortex l, if it were νε is given by 


Lm M: [2 n 


-u AV, ——- —— -—--—— |j 0) 
Z "Zr Z,-42 2-9 SUE 
| 2 z 


and likewise for vortex 2; 


fy = h IUE e (39) 
+27 "zz ni ni 


Next, the feeding aspect of the vortex cores is handled by consider- 


—U FAV, = = IN 








ing the vortex sheet, or feeding layer, as the combination of a single 
concentrated vortex, situated at the center of vorticity, with a con- 
necting sheet to the body of vanishingly small vorticity. This sheet 
can now be considered as a mathematical branch line between 2. апа δι 


for vortex 1, and Z, and 2, for vortex 2. Since there is a discontinu- 


b 
ity, conceptually, στους these branch lines, there is a velocity 


difference and therefore a pressure difference. This pressure difference 


14 
A free vortex is one that is no longer being fed circulation. 
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2 ιν 


Feeding sheet 1 


х 





Feeding sheet 2 


Figure !2. Model A, Indicating Capability for 
Asymmetry. 
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can be shown ([43] and [44]) to be equal εν р. . Since the dís- 


tances over which the pressure differences act are (2, -2,,) and(Z, -Z, 
for vortices 1 and 2 respectively, the forces acting perpendicular to 
these sheets are given by ipl" (2, -2,) and cp Α (2,-2,). The 


forces acting directly on the vortices is given, for vortex l, by 


, AWe) 2 
ipl, ( д 2. нь 2, ) ) 
Te 
and for vortex 2, by 


Í aa 7 


d W(2) | 


Af la( Ze EE 


2222 

Thus to make the net forces on each sheet and its associated concen- 
trated vortex πο... the sum of the respective forces must equal 
zero. Therefore, algebraically summing each pressure force with its 


proper vortex force, cancelling out i апа 2 , and rearranging the terms, 


the two expressions formed become; for vortex 1: 


e 











n d WiZ) 
Z,+(2,—Z )— = — | ) (40) 
uw dZ lj. 
and for vortex 2; 
W(2) en 


Г 
Z,t(Z, 2) = m 


These equations indicate that unless E is zero, the nth vortex does 
not move with the local fluid velocity, u + iv › but with an "attached" 


velocity of 208 


1 
>The combination of vortex sheet and the concentrated vortex in 


a fluid medium has to be force free since no external force is exerted 
on the system. 
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where the constants K are defined by: 
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from equation (29). 
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e > 
Since the M vortices are unattached, Z, and Z, are calculated 


Since there are more unknowns than there are equations, the system 
is still ill-defined. Four more equations can be gotten by imposing 
the Kutta condition at the point from which the vortex cores are fed. 
This is done by assuming that the fluid velocity relative to the cylin- 
der at these points is zero, i.e., Z, and Z, are stagnation points. Ву 


T b 


solving equation (36) for -u, t iv, and for -u, - iv, , then setting 


T T b b' 


these equations equal to zero and rearranging in matrix form, the two 


new equations become; 


/ АК | |; [к (К, + и) | 


\ 


27 | "m (42) 
\ 


-A, Ks / р | КК, + 27) 


Ej 

The Ks are the same ones associated with equation (43). Though not in- 
cluded in equations (38) and (39) for the sake of simplicity, the above 
equations include the extra summation terms, K), and Kig?’ necessary for 
the calculation if more than two vortex cores are present. To give the 
velocity equations (38) and (39) the same versatility, all that is neces- 
sary is to multiply these terms by 1/27, adding the first one to right 
hand side of (38) and the second one to (39). This has been done in 
the computer program for this model. 

The next two equations are gotten by solving the first equation of 


i 
(42 for Г and the second one for | . Then the time derivatives, [3 


2 
and j» are taken, assuming that τ and Z, are constants. Since the 
procedure is long and involved, but quite straight forward, only the re- 
sults will be presented; see equations (43). These equations also in- 
clude the terms necessary for more than two vortex cores. There is also 


a provision made in these equations for 02 to be taken into considera- 


tion. This quantity was set to zero for the study undertaken in this 
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thesis, but could be utilized in future work. Note that the equations 
e 

have been non-dimensionalized. This means that U. is also a non- 

dimensionalized quantity. 


-— 


i j Z, Z.. Z., Ž., and Z 
This now leaves three six unknowns, Z, Z. 1 2» 2) 2 
in the svstem of equations without a corresponding number of defining 


equations. The latter two unknowns can be easily handled by forming 


the conjugates of both sides of (40) and (41) to yield, 


— 5 A d wz) 





РА +(2,- 2 )— = $ 
l T [Г e 2-2, 
and , (44) 
z uz zh o dW 
2 тоер Г dz z=Z, 


as mentioned earlier, Z, and Z, are considered constants. This is an 
assumption based on experimental observations at higher laminar 

Reynolds numbers. The origins of these feeding sheets are in realitv 
functions of time, but theír motions are so small that the assumption 


seems justified. Besides, to consider Z, and 2 anything but zero, 


b 
would have made the derivation of (43) at least another order of magni- 
tude more tedious. This may, mods prove to be a necessary future 
sophistication of the model. As a result, these two locations are 
taken from experimental data to be about 50 degrees either side of the 


cylinder's downstream axis; i.e., O,, = 50 and e = -50 degrees. 


T 
No mention thus far has been made concerning the propagation as- 
pects of this solution, however since time dependent solutions generally 
imply that starting and ending conditions are necessary, it is well to 
discuss these now with respect to Σι and 2. With the assignment of 


any values to these two parameters, there have been enough equations 


derived and assumptions made to have the proper balance between unknowns 
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and defining equations to describe the system at any given instant. Not 
just any value of Ζι and 2, is desirable however, since this will not 
give even an approximate picture of core growth and drag buildup. Bryson 
[29], through certain mathematical considerations, found that the path 
line of the vortex cores, as they grow and move away from the cylinder, 
is unique in the sense that for any given feeding sheet origin, Zn or Z ° 
there can be only one path, per vortex, followed as time progresses. 
Without going into the details of his work, he also found that the feed- 
ing points, like the FOppl curve point sah are mathematical equilibrium 
positions for symmetric vortex pairs. On the basis of these findings, 
it is also determined that cores can leave the cylinder's surface only 
in a direction 30 degrees downstream from the tangent to the cylinder at 
the sheet origin. Bryson arbitrarily selected the starting values for 
Ζι and 2, as a distance of 0.05c away from the origins along these re- 
spective preferred lines of departure. This practice was continued in 
this study even though the model can handle asymmetry and more than two 
vortex cores. It was felt that these are good assumptions for an in- 
itial study of asymmetry. 

Except for determining when and how the vortex sheets are cut, there 
are no more assumptions to be made in completing the solution. The pro- 


cedure for solving the previously developed equations and the develop- 


ment of the κ, + iF) equation are not complicated and follow next. 


16 sty š ; : 
The Foppl curve is a symmetric curve about the downstream X-axis 


Originating at the rear stagnation point, on which a vortex pair can be 
symmetrically located such that they will maintain their positions. This 
is much the condition that exists at Reynolds numbers below about 50. 

See [6, p. 223]. 
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Equations (40), (41), (43), and (44) can be rearranged and com- 
bined in a straight forward manner for the matrix equation (45). Both 
this matrix equation and equation (42) are solved using the partitioned 
matrix technique as given in Fox [39, p. 121]. Briefly, this involves 
breaking up the basic matrix equation into its real and imaginary parts 
and forming a new equation whose matrices are twice the order of the 


original. To describe the technique, suppose the original equation is; 


| ij = ich. 
where [A] is a square matrix of order j, while $B} Р and {C} are column 


matrices with j elements each. Since all the elements of A, B, and C 


are assumed to be complex, the matrices can be written as; 


[A] |^] +1 [^] 
5B Ly Я $8, } 
fc} fc, ] οι je; 


Substituting these expressions back into the original equation, results 


in 


i ! 
NE + i = i 
Carrying out the indicated multiplication, and separating the real and 


ὅ τε 
ΠΠ P [ ag) 2514 2 C5 | 


These two equations are combined to form the new matrix equation: 


KS 


Л 
jer 


the imaginary parts of both sides, yields 


—[ ss 
> 
(< 
AL, 
ΚΠ 
oe 
t 
= 
> 
=! 
Ln 
C 
= 
| 
Мс 
i 


and, 


(46) 
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Since there are no oe numbers involved, this equation is solved 
by existing numerical methods, i.e., the inverse of the coefficient matrix 
premultiplies the constant matrix to yield the unknown values of fB? д 
In the computer program for this model, the inverse is found using a 
standard inversion subroutine with row pivoting. 

The propagation aspect of the solution is handled in the program 
through the use of a simple Runga-Kutta integration technique [40, p. 
232] by finding the subsequent vortex positions, after a /At time step, 
from equations (40) and (41). There are, of course, more sophisticated 
techniques, but in view of the parameters assumed, it is felt that the 
increased accuracy is not necessary at this point. In fact, there is 
not really that much difference in results if a straight Eulerian ap- 
proximation is used. 

The final relation to be developed is that for Ch + ic, - This is 
done by completing the derivation of equation (30), or more specifically 
by finding the appropriate solution for equation (33). For this case, 
in which the cores are fed by a single feeding sheet originating at a 
point, Sarpkaya, [33], and [35], has previously derived the necessary 
equation, so the details will not be given here. The result, in non- 


dimensional form and adapted to this model is: 


: : > š l - 
Gori = eN +i (z-z) az] 


Fla -Z) L(Z,72;2)] n 


M 
+4 el. It LV, (0; „+4)) + 
mz 
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The image velocities for the free and attached vortices are both calcu- 
lated, in a sense, using equation (37). However, where the results of 
equation (29) are used directly for the vortices, the results of solving 
matrix equation (45) must be used for the attached vortices. The modifi- 


cation to equation (37), used for this purpose is; for the image of 





vortex l, = 
м 
а 5 
= 
and for that of vortex 2, — (48) 


Ζ = — £ e 
zi — 2 
2 
The computer program written for this model is listed in its 





entirety in the last portion of this thesis. It is written in FORTRAN 
IV as applied to the IBM 360/67 digital computer at the Naval Post- 

graduate School. In general, the calculations are carried out using 
sixteen digit numbers. The typical run, describing the impulsively 

И 
С 


started flow up to 2215.0, requires less than three minutes of 





computer execution time. Provided that the initial starting points, 2, 


апа 2, are accurately defined, the solution appears to be stable and 
well behaved. The stability problem was previously discussed in con- 
nection with the vortex path line being unique. If the initial estimate 
for Zi and 22 are even slightly off, the solution technique attempts to 


return to the path, but as a consequence there may be some abnormally 
ә 9 

large absolute values for JN los: 2, апа Е To increase the 

stability, smaller time steps are taken at the start and just after a 

vortex has been shed and a new vortex starts to grow. The normal time 


step is 0.05, while this reduced time step is 0.02. The program will 


accept as many as twenty vortices before it terminates execution. 
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Since anyone desiring to utilize this program must become intimately 
familiar with its internal workings, it is felt that any further dis- 
cussion is unnecessary here. There is a brief, but more detailed des- 
cription in Appendix A. 

It was possible, through the use of this program, to generate a 
reasonably good drag curve using asymmetric considerations. To start 
with, if no shedding was allowed, it was found that by multiplying the 
value of 4 during each time step, from equation (42) by 1.05, prior 
to its use in equation (45), that the value of D never changed sign. 
It did however, peak at about the same Aut that it had previously 
(no asymmetry), but instead of decreasing to zero and subsequently chang- 
ing sign; it decreased somewhat and then increased beyond bound. The 
drag curve followed somewhat the same pattern. By arbitrarily cutting 
the vortex sheet at the minimum Г values, the drag curve in Figure 13 
was generated. The computer listing for Model A, though meant to be 
modified for other uses, is in fact, the same program that produced 
this curve. It is realized that the amount of asymmetry imposed is much 
too large, but is an illustration of the possibilities. A good begin- 
ning point for future work, would be with von Karman's stability condi- 


tions or to investigate the works done by Roshko [19] and Bearman [27]. 


De MODEL B 

Since this model is an attempt to approximate the vortex sheets 
generated by boundary layer separations, it is much more sensitive to 
the placement of the nascent vortices than is Model A to the location 
of the feeding sheet origins. Hence, the problem of mathematically 
describing vortex sheet formation is more obvious. This model has a 


distinct advantage over the other, in that it is not necessary to 
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incorporate any assumptions into the derivation of the basic equations 
of motion that take the place of matrix equation (45). The assumptions 
are present, however, in the determination of the strength and initial 
nascent vortex positions. Probably the most important single advantage 
of this model, is that in general, once the separation points and vortex 
strengths have been determined, the vortex sheets will roll up to form 
vortex cores, taking care of the vorticity transport and diffusion auto- 
matically. To discuss these aspects further will require the develop- 
ment of the applicable equations, so this follows next. 
Since the nascent vortex is assumed to gain its total strength 

during the time step in which it is born, equation (29) can be used to 


determine its convective velocity. It is repeated here for convenience. 





— +X V= И „(2-1)+ Sle > 2- - ὃν ας. 2-2 eee г) (29) Вереає 


The vortex convective velocity, as indicated in k. '. ἢ discussion 
of equation (29), is calculated by setting Z = Z, and ignoring the term 
in the first summation that has the denominator equal to zero. The 


velocity of the image vortex can still be calculated from 





i -Un tt Vy 
Uoan ha и = —5 e (37) Repeat 
=n 


The drag and lift coefficients are again calculated from equation (30), 
however, this time, some liberties can be taken with its use. Again, 
since the nascent vortex is born during one time step, a part of the 
effect of the unsteady portion of (30), more specifically parts of equa- 
tion (33), can be assumed negligible. In fact, a study of the results 
of [41] indicates that the algebraic sum of the values of the terms 


concerned with feeding and growth of the vortices actually only account 


63 


ἘΝ 


Birt^ Area 
| 







Birth Area 


T 
KU 







+H 


Birth Area 
Birth Area 4 
2 
rin E: Configuration of Model B Showing tne 


Birth Areas and the Rotation Direction 
of the Vortices Born Therein. 


0G 


for about three percent of the drag. The two terms being discussed are 
equivalent to the terms inside the first and second sets of brackets on 
the right hand side of equation (47). Since in this case M = N, the rest 


of the terms become applicable to this model. Therefore, 


| N 
CC, = CHU, ὍΣΩΝ VL, 3-4 V, —(U; +4 78) e (49) 
N=!) 

The system of basic equations is now complete, and like Model A, 
it has been programmed to allow asymmetry. The programming of equation 
(49) does not include the U. term since during the study conducted 
for this thesis, only a uniform, steady cross-flow of О ] was 
considered. It is simple enough to include this term, but then the 
deleted terms may become more important such that some kind of deter- 
mination might have to be made for them. It would be hoped that by 
the time these terms would be needed, the uniform flow problem would 
be solved and its results would be used to make this determination. 

This program is also written in FORTRAN IV, but unlike that of 
Model A, its calculation is carried out in only seven digits. The 
program, as written, will allow as many as one hundred nascent vortices 
to be born in any of the four birth areas during any given run. A run, 
which goes through the complete 100 time steps allowed, requires approxi- 
mately 30 minutes of computer execution time. Since the propagation 
aspects of the model are handled by Eulerian approximations, this could 
have serious error implications. However, since an approximate solution 
is all that is presently desired, it is felt that these errors are ac- 
ceptable. The large execution time is due primarily to the large number 


of vortices involved in the evaluation of (29), particularly at high Aet, 
τ 
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This must be done for each vortex, during each time step. The model also 
has the option of considering only the primary birth areas, 1 and 2, 
rather than all four. 

Before discussing some of the output from this model, the general 
method of determining vortex strengths must be considered. It has been 
shown by Fage and Johansen [42] that in the absence of back-flow and in 
a steady flow situation, the amount of vorticity leaving a separation 


zone can be accurately approximated by 


9[ ° 
πω --- le.sv'] ) (50) 
ens 

where V is the fluid velocity just outside the boundary layer at the 
point of separation. Thus, making the time steps, /At, small enough, 


the nascent vortex strength can be given by 
2 
\ = Jo. 5 V At] ә (51) 


The key problem now is to decide where and how to measure V and where 
the nascent vortex of this strength is placed. The four birth area 
model purports to take care of back-flow through the use of the second- 
ary birth areas, 3 and 4. Now since the maximum velocities on or near 
the cylinder's surface in these regions are very small as compared to 
the velocities in the primary regions, the effect of making some as- 
sumptions here is not too critical. The assumption made in most of the 
computer runs was that the maximum back-flow velocity on the cylinder's 
surface in that birth area was used in equation (51). From experimental 


studies, these back-flow vortices appear to separate at points located 
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Figure 15 
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Drag coefficient curve generated 
by Model B for an asymmetric case. 
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about + 55 degrees from the positive x-axis on the cylinder's surface; 
hence, these are the assumed angles for the model. Instead of placing 
them on the surface, however, the model was generally programmed to 
place them at about O.Olc. This still leaves the problem of placing 
the primary nascent vortices; to the author's knowledge, no truely 
satisfactory method does exist. The drag curves, and flow patterns 
shown in the figures of this subsection were all generated using arbi- 
trarily selected placement locations. The maximum velocity found on 
cylinder's surface in these birth areas was, however, used to compute 
the strength from (51). For all the birth areas, the sign given to the 
resultant value of strength is as depicted in Figure 14. This sign con- 
vention comes from boundary layer considerations and the direction of 
the maximum velocity in the respective birth area. 

In order to prevent the extremely high convective velocities that 
can occur if and when two vortices get too close together, the model 
provides for coalescence and cancellation of the vortices depending on 
their comparative strengths. The vortices are combined if the alge- 
braic sum is not zero, in which case they are cancelled, by creating a 
new vortex having that sum for its strength and located at the "center 
of gravity" of the two vortices. With a great many vortices this pro- 
cedure can consume much computer execution time. 

As in the case of Model A, any future user of this program would 
of necessity have to have an intimate knowledge of its inner workings. 
Appendix A contains a brief description of the parameters and the major 
sections of the program. 

This model was used to test a great many ideas using a great many 


runs. Figures 15 through 19 are presented to illustrate the kind and 
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type of output the model will generate. The program included with this 
thesis generated these figures. The primary vortices are injected into 
the stream at points whose angles from the positive x-axis are τ» 


degrees, and are located 0.045c from the cylinder's surface. 
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IV. DISCUSSION OF RESULTS AND CONCLUSIONS 


The primary result of this study is that reasonably representative 
drag coefficient curves can be generated using the two mathematical 
models developed. This is particularly true in the case of Model B, 
where the curves can be generated for both symmetric and asymmetric 
wake development. Figures 15 through 19 show the drag curve and wake 
development for an asymmetric wake case, while Figures 20 through 24 
show the symmetric wake case. The assumptions involved only the place- 
ment locations of the nascent vortices and the coalescence distance used 
to prevent abnormally high velocities if two nascent vortices come too 
close together. The symmetric wake case was generated by making the same 
assumptions for placement location as in the asymmetric case, however, 
the coalescence distance was increased. Where this distance was 0.06с 
for asymmetry, in the symmetric case it was 0.08c. These values were 
arbitrarily selected. 

In the case of Model B, the curve in Figure 15 was generated by 
assuming asymmetric generation of vorticity and assuming the vortex 


Г | reached its first post-peak 





sheet is cut, freeing the core, when 
minimum. The amount of asymmetry and the time of sheet cutting is 
loosely based on the results of Sarpkaya's experimental work [24]. 
Aside from the above results, numerous computer runs were made in 
an effort to develop a method by which the separation characteristics 
could be analytically described. Though these attempts were not en- 
tirely successful, some information was gained about the character of 
the vortex sheet origin and the requirements necessary to produce a 


reasonably realistic drag curve and wake structure. Of the two models, 
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Figure 20 
fficient curve generated 


y Model B for a symmetric case. 
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Figure 23 


Wake development corresponding 
to U,t - 7.4 in Figure 20. 
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Figure 22 


Wake development corresponding 
to Ut = 4.6 in Figure 20. 
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Model B was the most illuminating since the concept of a string of small 
strength vortices is more nearly like nature than is the concept upon 
which Model A is based. As a result, Model B was utilized more exten- 
sively in this study. It is for this reason then, that the majority of 
the remaining remarks in this section will be concerned with results from 
this model. 

Of the various factors that go into the final description of the 
wake and computation of a drag curve, it was found that the location at 
which the nascent vortices are placed in Model B is by far the most 
critical parameter involved. Of the two items which determine the loca- 
tion of the vortex, angle and distance from the cylinder's surface, it 
appears that the latter is more important. It was observed that there 
was a tendency for the primary vortices to 'hug' the surface (Figure 25), 
all the way to the rear stagnation point, if their initial placement 
distance was too small. "Тоо small' in this case appears to mean less 
than about 0.03c. In Figure 25, the nascent vortices were placed at 
0.01c away from the cylinder's surface, at angles of * 95 degrees. In- 
creasing or decreasing IF| did not seem to have much effect on this 
tendency. Changing the secondary vortex birth location, for a fixed 


' nor did it appreci- 


primary location, did not prevent this "hugging,' 
ably change the roll-up tendency of the primary vortices if they were 

being shed reasonably well. Though the resulting wake formation in 

Figure 25 is not that desired for Reynolds numbers less than about 2 x one 
there is a striking resemblance to real flow patterns evident in the 

super critical Reynolds number region. By making approximately 50 com- 
puter runs in which only the birth location of the primary vortices was 


varied, it was found that a distance of about 0.045c and an angle of 


about + 95 degrees from the rear stagnation point, produced the most 
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reasonable flow patterns and drag curve; Figures 15 through 19. This 
shedding angle is comparable with the observed angle of separation of 
about 98 degrees (measured in the same manner). The distance from the 
cylinder's surface however, is about twice the value predicted by the 
Blasius series solution of the boundary layer equations. It is suggested, 
therefore, that any future work in this area might well begin with a 
comprehensive study of the boundary layer transition region. А соп- 
sideration of these regions above flat plates could be the starting 

point since the added problem of body curvature is not present. 

It does not seem possible, after studying the works of many in- 
vestigators, that Thompson [13] could have produced asymmetry with 
equations that are in themselves symmetric, were it not for the very 
real problem of computer round-off error. The fact that Figure 15 ex- 
hibits asymmetry would be of deep concern had Abernathy and Kronauer 
[10] not found that strings of small strength vortices would in fact roll- 
up in an asymmetric manner. In the case of Model B, the predominate 
errors are due to assuming Eulerian propagation approximations and not 
from round-off. The Eulerian approximations are, of course, symmetric, 
so it is assumed that the asymmetry is in fact due to the interactions 
Abernathy and Kronauer observed. Interestingly enough, when using Model 
B asymmetry was not always produced, which leads to another conclusion; 
there is a critical combination of separation points and spacing dis- 
tance between sheets to cause the necessary interaction that forms the 
Karman street. On the other hand, reasonable drag coefficient curves 
could be generated without any asymmetry present. The overall conclusion 
is that the drag curve is predominatly influenced by the rate of vorti- 
city creation and its subsequent diffusion into the wake. This then, 


points right back to the over-riding problem of boundary layer separation. 
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Since Model A feeds large vortex cores rather than creating vortex 
sheets, it would appear this model is better suited for describing the 
flow at Reynolds numbers less than about 50. By increasing the angles 
of the sheet origin, the drag curve peak can be increased to almost any 
desired value. In order to make the model applicable to the higher 
Reynolds number range, it will be necessary to find some scientifically 
based method of shedding these cores. The stability criterion of von 


Karman and others may provide the answers. 
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APPENDIX A 


А. DESCRIPTION OF COMPUTER PROGRAM FOR MODEL A 

The complete program for this model includes a primary or main pro- 
gram with seven subroutines that accompany it. There is one other sub- 
routine that is called that is not part of the listing included with 
this thesis. The subroutine is PLOTP and is called from the resident 
subroutine library of the Naval Postgraduate School Computing Center. 
The subroutine is used to produce printer plots, such as Figures 15 
through 19, instead of using the CALCOMP continuous line plotters as 
output, so it is not critical to the actual program calculations. The 
curve plotting subroutines or procedures of other systems could be 
easily substituted for this one. 

At the beginning of the main program and each subroutine is a short 
series of comment cards that indicate the purpose of the program and 
generally the parameters used therein. In general, the same notation, 
in FORTRAN format, is used in the program as was used to develop the equa- 


tions in Section III(c); for example K, in equation (43) is coded as Kl, 


1 
ог 2, is coded as Z1B. Each computational section of the main program 
is preceded by comment cards indicating what computations that section 
is performing. 

Since there is no on-line function for extracting the imaginary 
part of a double precision complex number in the IBM version of FORTRAN 
IV, EQUIVALENCE statements had to be used to recover these values. 

It is assumed that the reader already has a working knowledge of 


FORTRAN coding techniques so that a definition list of the various para- 


meter names used in the program will be sufficient guide for future use, 
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The need for understanding the program is brought about by the fact that 
in order to study vortex sheet cutting procedures, the logic of the 
coding will have to be changed each time a different procedure is used. 
The section of this thesis entitled Computer Output, contains several 
pages of output that are to be used as trial us checks for a new pro- 
gram deck punched from the program listing. 

The type specification and precision of the following parameters is 
clearly indicated in the first few statements of the program, so there 
is no neea for repeating it here. All arrays are stored in vector form, 
i.e., the values of the array elements are Score one column following 
another, going from left to right in the array. 

If a parameter is not found in the following list (equation para- 
meters described above are not included) but found in the program, then 
it is a dummy parameter used for intermediate operations only and there- 


fore has no special significance as compared to those listed. 


В. DEFINITION OF PARAMETERS 

AA (144) Vector containing the partitioned matrix elements. 

AC (12) Vector containing the partitioned elements of the right 
hand side of equation (45). First six elements contain 


the real parts, the second six the imaginary parts. 


AG(12) Vector containing the partitioned elements of the un- 
knowns in equation (45). 


CDCL(400) Vector containing the lift and drag coefficient values 
for each time step. 


CV (20) Vector containing the free vortex velocities; corresponds 
to ZF(20) and GF(20). 

DG1(401) Contains Г for each time step. 

DG2 (401) Contains Ё for each time step. 
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DG 1MAX 
DG2MAX 
DG1MIN 
DG2MIN ; 


DT 
DUDT 
EG 
EGS 


GF (20) 


GG1 
GG2 


IT 


INNB 


INNT 


NCHDT 


NFV 


NFVT 


NIT 


NKER(2) 


NPL 


: : с TL 
Contains the values of the maximum and minimum f., and [ж 


for a run; used to determine when to cut the vortex 
sheets. 


Contains value of Zt; time step size. 
U, 

-12 
A tolerance value, usually 1.0x10 : 
EG squared. 


Vector containing the values of the strengths ( o 
of the free vortex cores; corresponds to the position 
vector ZF(20). 


Value of ls during a time step. 
Value of E during a time step. 


Counter value for the pass number through the program. 
ES ae Cs t 


If bottom vortex is shed during a time step, this para- 
meter equals 0, if vortex is not shed, it equals l. 

It is used to determine whether a new vortex should be 
started from bottom sheet. 


Has same indication as INNB, but applies to the top 
vortex and sheet. 


Indicates the pass at which the value of DT is changed 
to prevent instabilities in the solution. 


Counter used to indicate the number of free vortices 
existing. NFV = 20 


Contains the maximum value of NFV during any run. If 
NFVT >20, the DIMENSION statement must be changed as 
indicated in the program. 


Contains the maximum value of IT during any run. If 
NIT > 400, the DIMENSION statement must be changed as 
indicated in the program. 


Vector containing the matrix singularity determination 
from subroutine INVERT. NKER(1) applies to equation 
(42) while NKER(2) applies to equation (45). 


Indicates number of passes through program between wake 
development plots. 
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2 


NPRT 


PI 
PI2 
T(401) 


ZA(36) 


ZAA (^) 


ZC(6) 


ZCC(2) 


ZCP (100) 


ZF (20) 


ZG(6) - 


ZGG(2) 


Indicates number of passes through program between 
status printouts. 


E 
ay 2 


Contains value of total elapsed == for each time step. 





Vector containing the values of the elements in the co- 
efficient matrix of equation (45). 


Same as ZA(36), but applies to equation (42). 


Vector containing the elements of the column matrix on 
the right hand side of equation (45). 


Same as ZC(6), but applies to equation (42). 


Contains the circle circumference points and scaling 
points used for plotting purposes. 


Contains the free vortex positions during a time step; 
corresponds to the vector GF(20) on a one for one basis. 


Contains the values of the unknowns in equation (45). 


Same as ZG(6), but applies to equation (42). 
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APPENDIX B 


А. DESCRIPTION OF COMPUTER PROGRAM FOR MODEL B 

This program is much shorter from a coding standpoint than is Model 
A. There is only one subroutine besides the plotting routine, PLOTP, 
needed for operation of the program, and it is included in the listing. 
The general comments previously made for Model A, apply equally well 
here. The only major difference, besides the concepts upon which the 
two are based, is that this program is written using single precision 


arithmetic rather than double precision as in the case of Model A. 


B, DEFINITION OF PARAMETERS 


AVN(4) Values, in degrees, of the nascent vortex placement 
angles actually used in the calculations. AVN(1) 
applies to birth area 1; AVN(2) applies to birth area 


εις. 
CDCL(100) Same as Model A. 
CLOCKR These parameters are used to determine the actual 
CLOCKS elapsed computer time for each pass and run; 
CLOCKT can be eliminated from the program without hurting 


its purpose. 


CVV (400) Contains the value velocities for each vortex during 
a time step. CVV(1)to CVV(100) Birth Area 1 
CVV(101) to CVV(200) Birth Area 2 
CVV(201) to CVV(300) Birth Area 3 
CVV(301) to CVV(400) Birth Area 4 


DISP Same as DT in Model A 

EG Same as Model A 

EGS Same as Model A 

GK(400) Contains the value of the vortex strengths during a 
time step; corresponds to CVV(400) on a one for one 
basis. 

IT Same as Model A 
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ITIME(0) 


MPA 


МР(4) 


М2р 


МВ 


МСАМ 


NCL 


NIT 


NP1 


NVM(4) 


PI 
RE 
TP(101) 


VELM(4) 


ZCC(100) 


An IBM on-line function used to make real time 
available to a FORTRAN IV program; can be eliminated 
without hurting the purpose of the program. 


The value of nascent vortex placement distance as a 
function of RE as generated by formulas from Sarpkaya's 
paper [35]. Used for comparative purposes. 


Contains the actual nascent vortex placement distance 
from the cylinder's surface used in the calculations. 
MP(1) corresponds to AVN(1) or NVM(1); MP(2) corres- 
ponds to AVN(2) or NVM(2); etc. 


Radius of the arc along which the search for maximum 
velocities is made; usually equal to 1.0. 


Minimum distance two vortices can be from each other 
without being coalesced or cancelled. 


Total number of birth areas; usually equal to 4, but 
if only the primary nascent vortices are being studied, 
then NB=2. 


Counter indicating number of vortices cancelled during 
a run. 


Counter indicating number of vortices coalesced during 
a run. 


Same as Model A. 


Indicates number of passes, times ten, between wake 
development plots. 


Counter indicating number of primary vortices generated 
in birth area 1 during a run. It is also equal to the 
number generated in birth area 2. 

Contains the values, in degrees, of the maximum velocity 


locations resulting from a search for same in the birth 
areas. NVM(1) corresponds to AVN(1); etc. 


7 

Reynolds number; used in conjunction with MPA. 

Same as T(401) in Model A. 

Contains the values of maximum velocity found during the 
search for same in the birth areas. VELM(1) corresponds 


to NVM(1); etc. 


Same as ZCP(100) in Model A. 


89 


ZCFP(51) Contains the cylinder's circumferential locations from 
135 to 85 degrees (measured from the rear stagnation 
point). 


ZCRP(61) Contains the cylinder's circumferential locations from 
O to 60 degrees (measured from the rear stagnation point). 


ZF (400) Contains the locations of the vortices during a time step. 
Corresponds to CVV(400) and GK(400). 
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